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Modeling  primary  breakup:  A three-dimensional 
Eulerian  level  set /vortex  sheet  method  for 
two-phase  interface  dynamics 

By  M.  Herrmann 


1.  Motivation  and  objectives 

Atomization  processes  play  an  important  role  in  a wide  variety  of  technical  applications 
and  natural  phenomena,  ranging  from  inkjet  printers,  gas  turbines,  direct  injection  IC- 
engines,  and  cryogenic  rocket  engines  to  ocean  wave  breaking  and  hydrothermal  features. 
The  atomization  process  of  liquid  jets  and  sheets  is  usually  divided  into  two  consecutive 
steps:  the  primary  and  the  secondary  breakup.  During  primary  breakup,  the  liquid  jet  or 
sheet  exhibits  large  scale  coherent  structures  that  interact  with  the  gas-phase  and  break 
up  into  both  large  and  small  scale  drops.  During  secondary  breakup,  these  drops  break 
up  into  ever  smaller  drops  that  finally  may  evaporate. 

Usually,  the  atomization  process  occurs  in  a turbulent  environment,  involving  a wide 
range  of  time  and  length  scales.  Given  today’s  computational  resources,  the  direct  nu- 
merical simulation  (DNS)  of  the  turbulent  breakup  process  as  a whole,  resolving  all 
physical  processes,  is  impossible,  except  for  some  very  simple  configurations.  Instead, 
models  describing  the  physics  of  the  atomization  process  have  to  be  employed. 

Various  models  have  already  been  developed  for  the  secondary  breakup  process.  There, 
it  can  be  assumed  that  the  characteristic  length  scale  l of  the  drops  is  much  smaller  than 
the  available  grid  resolution  Ax  and  that  the  liquid  volume  fraction  in  each  grid  cell  0; 
is  small,  see  Fig.  1.  Furthermore,  assuming  simple  geometrical  shapes  of  the  individual 
drops,  like  spheres  or  ellipsoids,  the  interaction  between  these  drops  and  the  surrounding 
fluid  can  be  taken  into  account.  Statistical  models  describing  the  secondary  breakup 
process  in  turbulent  environments  can  thus  be  derived  (O’Rourke  1981;  O’Rourke  & 
Amsden  1987;  Reitz  1987;  Reitz  & Diwakar  1987;  Tanner  1997). 

However,  the  above  assumptions  do  not  hold  true  for  the  primary  breakup  process. 
Here,  the  turbulent  liquid  fluid  interacts  with  the  surrounding  turbulent  gas-phase  on 
scales  larger  than  Ax,  resulting  in  highly  complex  interface  dynamics  and  individual  grid 
cells  that  can  be  fully  immersed  in  the  liquid  phase,  compare  Fig.  1.  An  explicit  treatment 
of  the  phase  interface  and  its  dynamics  is  therefore  required.  To  this  end,  we  propose  to 
follow  in  essence  a Large  Eddy  Simulation  (LES)  type  approach:  all  interface  dynamics 
and  physical  processes  occurring  on  scales  larger  than  the  available  grid  resolution  Ax 
shall  be  fully  resolved  and  all  dynamics  and  processes  occurring  on  subgrid  scales  shall 
be  modeled.  The  resulting  approach  is  called  Large  Surface  Structure  (LSS)  model. 

In  order  to  develop  such  a LSS  model  for  the  turbulent  primary  breakup  process, 
one  potential  approach  is  to  start  off  from  a fully  resolved  description  of  the  interface 
dynamics  using  the  Navier-Stokes  equations  and  include  an  additional  source  term  in  the 
momentum  equation  due  to  surface  tension  forces  (Brackbill  et  al.  1992).  In  order  to  track 
the  location,  motion,  and  topology  of  the  phase  interface,  the  Navier-Stokes  equations  are 
then  coupled  to  one  of  various  possible  tracking  methods,  for  example  marker  particles 
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Figure  1.  Breakup  of  a liquid  jet. 


(Brackbill  et  al.  1988;  Rider  & Kothe  1995;  Unverdi  & Tryggvason  1992),  the  Volume-of- 
Fluid  method  (Noh  & Woodward  1976;  Kothe  & Rider  1994;  Gueyffier  et  al.  1999),  or  the 
level  set  method  (Osher  & Sethian  1988;  Sussman  et  al.  1994,  1998).  Then,  introducing 
ensemble  averaging  or  spatial  filtering  results  in  unclosed  terms  that  require  modeling 
(Brocchini  & Peregrine  2001  a, b).  Unfortunately,  the  derivation  of  such  closure  models  is 
not  straightforward  and,  hence,  has  not  been  achieved  yet.  This  is  in  part  due  to  the  fact 
that,  with  the  exception  of  the  surface  tension  term,  all  other  physical  processes  occurring 
at  the  phase  interface  itself,  like  for  example  stretching,  are  not  described  by  explicit 
source  terms.  Instead,  they  are  hidden  within  the  interdependence  between  the  Navier- 
Stokes  equations  and  the  respective  interface  tracking  equation.  Thus,  a formulation 
containing  the  source  terms  explicitly  could  greatly  facilitate  any  attempt  to  derive  the 
appropriate  closure  models. 

To  this  end,  a novel  three-dimensional  Eulerian  level  set/vortex  sheet  method  is  pro- 
posed. Its  advantage  is  the  fact  that  it  contains  explicit  source  terms  for  each  individual 
physical  process  that  occurs  at  the  phase  interface.  It  thus  constitutes  a promising  frame- 
work for  the  derivation  of  the  LSS  subgrid  closure  models. 

This  paper  is  divided  into  four  parts.  First,  the  level  set/vortex  sheet  method  for  three- 
dimensional  two-phase  interface  dynamics  is  presented.  Second,  the  LSS  model  for  the 
primary  breakup  of  turbulent  liquid  jets  and  sheets  is  outlined  and  all  terms  requiring 
subgrid  modeling  are  identified.  Then,  preliminary  three-dimensional  results  of  the  level 
set /vortex  sheet  method  are  presented  and  discussed.  Finally,  conclusions  are  drawn  and 
an  outlook  to  future  work  is  given. 


2.  The  level  set/vortex  sheet  method 

The  aim  of  the  level  set/ vortex  sheet  method  is  to  describe  the  dynamics  of  the  phase 
interface  P between  two  inviscid,  incompressible  fluids  1 and  2,  as  shown  in  Fig.  2.  In  this 
case,  the  velocity  u,  on  either  side  i of  the  interface  F is  determined  by  the  incompressible 
Euler  equations,  given  here  in  dimensionless  form, 

V • ui  = 0 , 

^ + («1.V)ui  = -ivp, 

ot  pi 


(2.1) 

(2.2) 
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Figure  2.  Phase  interface  definition. 


subjected  to  the  boundary  conditions  at  the  interface  T, 


[(ui  ~ «2)  ’ n] 

[n  x (tt2  - ui)] 
\P2-Pl) 


= 0 
r 


r We* 


(2.3) 

(2.4) 

(2.5) 


and  at  infinity, 


lim  Ui  = ±«oo  . (2.6) 

v-’ioo 

Here,  n is  the  interface  normal  vector,  7)  is  the  vortex  sheet  strength,  and  k is  the  local 
curvature  of  T.  The  Weber  number  is  defined  as 


We  — pief%iief /YjLxe{ , (2.7) 

where  E is  the  surface  tension  coefficient  and  pre f,  uref , and  Lref  are  the  reference  den- 
sity, velocity,  and  length,  respectively.  An  interface  subjected  to  the  above  boundary 
conditions  is  called  a vortex  sheet  (Saffman  & Baker  1979). 

The  partial  differential  equation  describing  the  evolution  of  the  vortex  sheet  strength  77 
can  be  derived  by  combining  the  Euler  equations,  Eqs.  (2.1)  and  (2.2),  with  the  boundary 
conditions  at  the  interface,  Eqs.  (2.3)-(2.5),  resulting  in  (Pozrikidis  2000) 

di] 

— + u ■ V77  = — n x [(77  x n)  ■ Vu]  + n [(Vw  • n)  • 77] 

(71  x Vk)  + 2Ati  x a.  (2.8) 

Here,  A = (p\  — p2)/{pi  + P2)  is  the  Atwood  number  and  a is  the  average  acceleration 
of  fluid  1 and  fluid  2 at  the  interface.  The  major  advantage  of  Eq.  (2.8),  as  compared  to 
a formulation  based  on  the  Euler  equations,  is  the  fact  that  Eq.  (2.8)  contains  explicit 
local  individual  source  terms  on  the  right-hand  side  describing  the  physical  processes  at 
the  interface.  These  are,  from  left  to  right,  two  stretching  terms,  a surface  tension  term 
To-,  and  a density  difference  term. 

In  addition  to  the  evolution  of  the  local  vortex  sheet  strength,  Eq.  (2.8),  the  location 
and  motion  of  the  phase  interface  itself  has  to  be  known.  To  this  end,  vortex  sheets  are 
typically  solved  by  a boundary  integral  method  within  a Lagrangian  framework  where 
the  phase  interface  is  tracked  by  marker  particles  (Baker  et  al.  1982;  Pullin  1982;  Hou 
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et  al.  1997,  2001;  Rangel  & Sirignano  1988).  Marker  particles  allow  for  highly  accurate 
tracking  of  the  phase  interface  motion  in  a DNS.  However,  the  introduction  of  ensemble 
averaging  and  spatial  filtering  of  the  interface  topology  is  not  straightforward  and  hence 
a strategy  for  the  derivation  of  appropriate  LSS  subgrid  closure  models  is  not  directly 
apparent. 

Level  sets,  on  the  other  hand,  have  been  successfully  applied  to  the  derivation  of  closure 
models  in  the  field  of  premixed  turbulent  combustion  (Peters  1999,  2000).  Thus,  instead 
of  using  marker  particles  to  describe  the  location  and  motion  of  the  phase  interface,  here, 
the  interface  is  represented  by  an  iso-surface  of  the  level  set  scalar  field  G(x,t ),  as  shown 
in  Fig.  2.  Setting 

G(x,  t)  |r  = G0  — const , (2.9) 

G(x,t)  > Go  in  fluid  1,  and  G(x,t)  < Go  in  fluid  2,  an  evolution  equation  for  the  scalar 
G can  be  derived  by  simply  differentiating  Eq.  (2.9)  with  respect  to  time, 


dG 

dt 


+ it-VG  = 0. 


(2.10) 


This  equation  is  called  the  level  set  equation  (Osher  & Sethian  1988).  Using  the  level  set 
scalar,  geometrical  properties  of  the  interface,  like  its  normal  vector  and  curvature,  can 
be  easily  expressed  as 


n 


VG 

|VG| 


V • n . 


(2.11) 


Strictly  speaking,  Eqs.  (2.8)  and  (2.10)  are  valid  only  at  the  location  of  the  inter- 
face itself.  However,  to  facilitate  the  numerical  solution  of  both  equations  in  the  whole 
computational  domain,  r]  is  set  constant  in  the  interface  normal  direction, 


Vtj  • VG  = 0 , 


(2.12) 


and  G is  chosen  to  be  a distance  function  away  from  the  interface, 

Ho*,.-1-  <2-13> 

Equations  (2.8)  and  (2.10)  are  coupled  by  the  self-induced  velocity  u of  the  vortex 
sheet.  To  calculate  u,  the  vector  potential  np  is  introduced, 


Ar/>  = u> . 


(2.14) 


Here,  the  vorticity  vector  u>  is  calculated  following  a vortex-in- cell  type  approach  (Chris- 
tiansen 1973;  Cottet  & Koumoutsakos  2000) 

u>(x)  = [ r](x')6(x  - x')5  (G(x')  - G0)  \\7G(x')\dx' , (2.15) 

Jv 

where  6 is  the  delta-function.  Then,  u can  be  calculated  from 


u(x)  = [ 5(x  — x')  (V  x ip)  dx' . (2.16) 

Jv 

In  summary,  Eqs.  (2.8),  (2.10),  and  (2.14)  - (2.16)  describe  the  three-dimensional  two- 
phase  interface  dynamics  and  constitute  the  level  set/vortex  sheet  method. 


2.1.  Numerical  methods 

Numerically,  Eqs.  (2.8)  and  (2.10)  are  solved  in  a narrow  band  (Peng  et  al.  1999)  by 
a third-order  WENO  scheme  (Jiang  & Peng  2000)  using  a third-order  TVD  Runge- 
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Figure  3.  LSS  interface  definitions. 


Kutta  time  discretization  (Shu  & Osher  1989).  The  redistribution  of  rj  (2.12)  is  solved 
by  a Fast  Marching  Method  (Sethian  1996;  Adalsteinsson  & Sethian  1999),  whereas  the 
reinitialization  of  G (2.13)  is  solved  by  an  iterative  procedure  (Sussman  et  al.  1994; 
Peng  et  al.  1999).  The  interested  reader  is  referred  to  Herrmann  (2002)  for  a detailed 
description  of  the  numerical  methods  employed  in  the  level  set /vortex  sheet  method. 


3.  The  LSS  model  for  turbulent  primary  breakup 

The  basic  idea  of  the  LSS  model  is  to  split  the  treatment  of  the  primary  breakup 
process  into  two  parts.  All  phase  interface  dynamics  occurring  on  scales  larger  than 
the  local  grid  size  are  explicitly  resolved  and  tracked  by  a level  set  approach,  whereas 
interface  dynamics  occurring  on  subgrid  scales  are  described  by  an  appropriate  subgrid 
model.  Furthermore,  the  LSS  subgrid  model  has  to  separate  out  all  broken  off  subgrid 
scale  liquid  drops  and  transfer  them  to  a secondary  breakup  model. 

The  level  set  equation  describing  the  interface  location  and  motion  on  the  resolved 
scales  can  be  derived  by  first  introducing  appropriate  interface  based  filters  (Oberlack 
et  al.  2001)  into  the  level  set  equation  (2.10),  see  Fig.  3, 

dG 

— + S-VG  = 0.  (3.1) 

Here,  T denotes  quantities  on  the  resolved  (filter)  scale.  Furthermore,  the  mass  transfer 
rate  rhp  into  the  secondary  breakup  model  has  to  be  taken  into  account, 


''Ad  f^x 

rhP  = piSpAdg  = -npi— P(D)D3dD,  (3.2) 

where  sp  is  the  subgrid  primary  breakup  velocity,  Aq0  is  the  local  surface  area  of  the 
resolved  interface,  and  P(D ) is  the  droplet  diameter  number  distribution, 


(3.3) 


where  N is  the  total  number  of  drops.  Then,  the  resolved  scale  level  set  equation  reads 
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where  n is  the  normal  vector  of  the  resolved  scale  interface, 


~ VG 
n = — — . 
|VG| 


(3.5) 


To  describe  the  phase  interface  dynamics  on  the  resolved  scale,  their  effect  on  the  flow 
field  has  to  be  taken  into  account  by  an  additional  source  term  T in  the  momentum 
equation, 

T = E kS(G  - G0)n  + TSG s • (3.6) 

Here,  the  first  term  on  the  right-hand  side  describes  the  effect  of  surface  tension  forces 
due  to  the  local  curvature  k of  the  resolved  scale  interface,  whereas  the  second  term, 
T sgSi  accounts  for  the  effect  of  the  subgrid  scale  surface  tension  forces  on  the  resolved 
scale  flow  field. 

Thus,  the  yet  unclosed  subgrid  terms  of  the  LSS  model  requiring  modeling  are  the 
subgrid  primary  breakup  velocity  sp,  the  droplet  diameter  number  distribution  P{D), 
and  the  subgrid  scale  surface  tension  effect  Tsgs  • As  previously  indicated,  these  subgrid 
terms  are  to  be  derived  from  the  level  set/vortex  sheet  method.  Performing  DNS  of 
the  primary  breakup  of  liquid  surfaces  and  sheets  in  turbulent  environments  will  help 
to  identify  characteristic  regimes  of  the  turbulent  primary  breakup  and  their  dominant 
physical  processes.  These  can  then  be  quantified  using  the  explicit  source  terms  in  the 
77-equation  (2.8),  thus  providing  guidelines  for  the  derivation  of  appropriate  LSS  subgrid 
models. 


4.  Results 

In  order  to  both  validate  the  three-dimensional  level  set /vortex  sheet  method  and  to 
demonstrate  its  ability  to  perform  DNS  of  the  primary  breakup  process,  the  results  of  two 
different  cases  are  presented.  First,  the  calculated  oscillation  periods  of  liquid  columns 
and  spheres  are  compared  to  theoretical  results.  Then,  the  breakups  of  a randomly  per- 
turbed liquid  surface  and  sheet  are  presented. 

4.1.  Oscillating  columns  and  spheres 

To  validate  the  proposed  level  set/vortex  sheet  method,  the  calculated  oscillation  periods 
T of  liquid  columns  and  spheres  of  mean  radius  R = 0.25,  center  xc  = (0.5, 0.5, 0.5), 
amplitude  e = 0.05 R,  and  Atwood  number  A = 0 are  compared  to  theoretical  results 
(Lamb  1945).  The  initial  vortex  sheet  strength  in  both  cases  is  set  to 

r](x,t  — 0)  = 0.  (4.1) 

All  calculations  are  performed  in  a unit  sized  box  resolved  by  an  equidistant  cartesian 
grid  of  128  x 128  and  128  x 128  x 128  nodes,  respectively. 

Figure  4 shows  the  distribution  of  the  surface  tension  term  Ta  of  the  77-equation  (2.8) 
in  the  x-,  y-,  and  2-direction, 

Ta  = (n  x Vk)  , (4.2) 

for  the  oscillating  sphere  of  mode  number  n = 5 and  Weber  number  We  = 10  calculated 
at  t = 0.  As  the  shape  of  the  sphere  indicates,  TV  in  the  x-direction  is  a factor  of  roughly 
four  higher  than  Ta  in  the  other  two  directions,  leading  to  the  predominant  oscillation 
in  the  y- 2-plane. 
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Figure  4.  Distribution  of  the  surface  tension  term  TV  in  the  z-direction  (top),  y-direction 
(left),  and  z-direction  (right)  for  the  mode  n — 5 oscillating  sphere  at  t = 0 and  We  = 10. 


Figure  5 depicts  the  comparison  of  the  oscillation  period  for  the  oscillating  columns 
on  the  left-hand  side  and  the  oscillating  spheres  on  the  right-hand  side  for  two  different 
Weber  numbers.  As  can  be  clearly  seen,  agreement  between  simulation  and  theory  is  very 
good. 


4.2.  Liquid  surface  and  sheet  breakup 

To  demonstrate  the  capability  of  the  proposed  level  set/vortex  sheet  method  to  simulate 
the  primary  breakup  process,  the  temporal  evolution  of  both  a randomly  perturbed 
liquid  surface  and  sheet  are  simulated.  In  the  case  of  the  liquid  surface,  the  on  average 
flat  interface  located  at  z = 0 is  perturbed  in  the  z-direction  by  a Fourier  series  of  64 
sinusoidal  waves  in  both  the  x-  and  y-direction  with  random  amplitude  0 < e < 0.01  and 
random  phase  shift.  In  the  case  of  the  liquid  sheet,  the  two  on  average  flat  interfaces  are 
located  at  z = — B/2  and  z = +B/2  and  are  again  perturbed  by  two  Fourier  series  of  64 
sinusoidal  waves.  The  thickness  of  the  liquid  sheet  is  set  to  B = 0.1. 

The  initial  vortex  sheet  strength  for  the  liquid  surface  is  set  to 


and  to 


r)(x,t  = 0)  - (-1,0,0) 


r](x,t 


(-1,0,0) 
( 1,0,0) 


z > 0 
z <0  , 


(4.3) 


(4.4) 
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n n 

FIGURE  5.  Oscillation  period  T of  liquid  columns  (left)  and  spheres  (right)  as  a function  of  mode 
number  n for  varying  Weber  numbers  We.  Lines  denote  theoretical  and  symbols  computational 
results. 


in  the  liquid  sheet  case.  Both,  the  surface  and  the  sheet  simulation  were  performed  in  a 
x-  and  ^-direction  periodic  box  of  size  [0, 1]  x [0, 1]  x [—1, 1]  resolved  by  a cartesian  grid 
of  64  x 64  x 128  equidistant  nodes.  In  both  simulations,  the  Atwood  number  is  A = 0. 
The  Weber  number  in  the  surface  simulation  is  We  = 500  and  the  Weber  number  in  the 
sheet  simulation  based  on  the  sheet  thickness  is  WeB  = 100. 

As  depicted  in  Fig.  6,  the  surface  shows  an  initial  growth  of  two-dimensional  Kelvin- 
Helmholtz  instabilities  (t  = 1).  These  continue  to  grow  (t  = 3)  and  form  three-dimensional 
structures  ( t — 5)  resulting  in  elongated  fingers  ( t = 6.5)  that  finally  initiate  breakup 
(t  = 8.0). 

The  liquid  sheet,  depicted  in  Fig.  7,  also  exhibits  the  initial  formation  of  two-dimensional 
Kelvin-Helmholtz  instabilities  (t  = 1)  that  continue  to  grow  (t  = 3)  until  the  liquid  film 
gets  too  thin  and  ruptures  (t  — 5).  Individual  fingers  are  formed  that  extend  mostly  in 
the  transverse  direction  ( t = 8)  and  continue  to  break  up  into  individual  drops  of  varying 
sizes  (f  = 12). 


5.  Conclusions  and  future  work 

A Eulerian  level  set/vortex  sheet  method  has  been  presented  that  allows  for  the  three- 
dimensional  calculation  of  the  phase  interface  dynamics  between  two  inviscid  and  in- 
compressible fluids.  Results  obtained  with  the  proposed  method  for  oscillating  columns 
and  spheres  show  very  good  agreement  with  theoretical  predictions.  Furthermore,  the 
applicability  of  the  method  to  the  primary  breakup  process  has  been  demonstrated  by 
simulations  of  the  breakup  of  both  a liquid  surface  and  a liquid  sheet. 

In  addition,  the  LSS  model  for  turbulent  primary  breakup  has  been  outlined,  and  all 
terms  requiring  subgrid  modeling  have  been  identified.  The  proposed  level  set /vortex 
sheet  method  has  the  advantage  that  it  allows  for  the  detailed  study  of  each  individual 
physical  process  occurring  at  the  phase  interface.  It  thus  provides  a promising  framework 
for  the  derivation  of  the  LSS  subgrid  models. 

Future  work  will  focus  on  including  the  effect  of  non-zero  Atwood  numbers  and  on 
coupling  of  the  level  set /vortex  sheet  method  to  an  outside  turbulent  flow  field.  Also,  the 
level  set /vortex  sheet  method  will  be  parallelized  making  use  of  the  new  domain  decom- 
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Figure  6.  Temporal  evolution  of  the  three-dimensional  liquid  surface  breakup,  A = 0, 

We  = 500. 


position  parallelization  of  the  Fast  Marching  Method  presented  in  Herrmann  (2003).  This 
will  allow  for  efficient  DNS  of  the  primary  breakup  process  to  help  identify  the  different 
regimes  of  turbulent  primary  breakup  and  their  dominant  physical  processes,  facilitating 
the  derivation  of  the  LSS  subgrid  models.  Finally,  combining  the  LSS  model  to  spray 
models  describing  the  secondary  breakup  will  allow  for  the  first  LES  of  the  turbulent 
atomization  process  as  a whole. 
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Figure  7.  Temporal  evolution  of  the  three-dimensional  liquid  sheet  breakup,  = 0, 

Wefl  = 100. 
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